Nonadiabatic dynamics of a Bose-Einstein condensate in an optical lattice 
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We study the nonequilibrium dynamics of a Bose-Einstein condensate that is split in a harmonic 
trap by turning up a periodic optical lattice potential. We evaluate the dynamical evolution of 
the phase coherence along the lattice and the number fluctuations in individual lattice sites within 
the stochastic truncated Wigner approximation when several atoms occupy each site. We show 
that the saturation of the number squeezing at high lattice strengths, which was observed in recent 
experiments by Orzel et. al, can be explained by the nonadiabaticity of the splitting. 
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Ultra-cold atomic gases in periodic optical lattice po- 
tentials have recently attracted considerable interest and 
inspired experiments, ej* on the Bose-Einstein conden- 
sate (BEC) coherence 111 |2( , superfluid dynamics 
the number-squeezed |(| and the Mott insulator (MI) 
states, and quantum information applications 0,0. An 
optical lattice provides a clean many-particle system with 
enhanced interactions, resulting in a unique opportunity 
to study strong quantum fluctuations. While the classical 
mean-field theories, such as the Gross-Pitaevskii equation 
(GPE), have been successful in describing the full multi- 
mode dynamics of weakly-interacting BECs, they have 
severe limitations in optical lattices, as they disregard 
thermal and quantum fluctuations, decoherence, and the 
information about quantum statistics. In this paper we 
study matter wave dynamics beyond the GPE by con- 
sidering a harmonically trapped finite-temperature BEC 
that is dynamically split by an optical lattice potential. 
We show that the experimentally observed saturation of 
the number squeezing at high lattice strengths [(| can be 
explained by the nonadiabaticity of the loading of atoms 
into the lattice. The thermal and quantum fluctuations 
are included within the truncated Wigner approximation 
(TWA). The multi-mode TWA provides a natural rep- 
resentation for the dynamical fragmentation of the ini- 
tially uniform BEC in the lattice and the transition to the 
regime that can also be described by the discrete Bose- 
Hubbard Hamiltonian (BHH). Moreover, the resulting 
highly occupied number squeezed states are also of grea t 
interest in the Heisenberg limited interferometry • 

We study the loading of atoms into the lattice within 
the TWA. The TWA may be obtained by using the 
familiar techniques of quantum optics [ill to de- 
rive a generalized Fokker-Planck equation (FPE) for the 
Wigner distribution of the trapped multi-mode BEC E3 . 
The TWA consists of neglecting the dynamical quantum 
noise, acting via third-order derivatives in the FPE, and 
results in a deterministic equation for the classical field 
4>w which coincides with the GPE: 



where C = — fi 2 V 2 /(2m) + V. The thermal and quan- 
tum fluctuations are included in the initial state of ipw 
in Eq. Q which represents an ensemble of Wigner dis- 
tributed wave functions. The neglected terms are small 
when the amplitudes of the Wigner distribution are 
large. The TWA and closely related approaches have 
previously been successful in describing atomic BECs 
El El El El IH E3 and optical squeezing 0. In 
particular, the TWA is shown to produce correctly, e.g., 
the Beliaev-Landau damping [15j and it has been argued 
that the TWA more generally provides an accurate de- 
scription for the short-time asymptotic behavior of the 
full quantum dynamics E3- 

We consider a BEC in a tight elongated cigar-shaped 
trap, with the trap frequencies uj = uj x <C uJ VtZ = uj±, 
and ignore the density fluctuations along the transverse 
directions. This results in an effective ID GPE for 
ipw(x,t) with g — (?id = 2thu±a in Eq. iJTjl. where a 
denotes the scattering length. The BEC is initially as- 
sumed to be in thermal equilibrium in a harmonic trap 
Vh(x) — mu) 2 x 2 /2. A self-consistent calculation of the 
initial state would involve solving the coupled Hartree- 
Fock-Bogoliubov equations for the condensate and non- 
condensate populations |2fj| . Here we resort to a simpler 
Bogoliubov approximation and expand the field operator 
^{x, t = 0) in terms of the BEC ground state amplitude 
ctQipo, with (djao) = Nq, and the excited states: 



j>0 



(2) 



where Uj(x) and Vj(x) (j > 0) are obtained from 



(C- n+ 2iV 3ii)|V'o| 2 ) Uj - N g 1D t/j 2 Vj = ejUj, 

*2„ 



(£ - n + 2N o9lD \vb \ 2 ) V j - NogmVo 



ejVj. (3) 



id t il>w = £il>w + g\ipw\ 2, >Pw , 



(1) 



Here &j are the quasiparticle annihilation operators, with 
(OjCtj) = fbj = [exp(/3ej) - l] -1 , (3 = l/k B T, and ip 
is ground state solution of the GPE with the chemical 
potential /i. 
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In the Wigner description we replace the quantum op- 
erators (ctj, &j) (for j > 0) by the complex random vari- 
ables (ctj,a*), obtained by sampling the corresponding 
Wigner distribution of ideal harmonic oscillators in a 
thermal bath 11]: 

W(aj, a*,) = - tanh (6) exp \-2\aA 2 tanh (6)1 , (4) 

7T 

where £j = /3ej/2. The Wigner function is Gaussian dis- 
tributed with the width fij + \. The nonvanishing contri- 
bution to the width at T — for each mode represents the 
quantum noise. The Wigner function returns symmetri- 
cally ordered expectation values, so (a*aj)w = % + 
and {otj)w = (a*)w — {ot 2 )w = 0, etc. 

For large BECs, iVo ^ 1, the main contribution to 
the matter wave coherence in the superfiuid regime is 
due to the thermal and quantum fluctuations of low- 
energy phonons and the quantum fluctuations of the ini- 
tial harmonically trapped BEC mode are not very im- 
portant. Consequently, we could treat the BEC mode 
do even classically. However, here we assume it to be 
in a coherent state and sample the quantum fluctua- 
tions according to the corresponding Wigner distribu- 
tion [H|: W(a ,a^) = 2exp(-2|a - W 1/2 | 2 )/7r, so that 

(ao)w — Nq^ 2 and (agao)w = No+h- Since we compare 
the matter wave coherence between the atoms in differ- 
ent lattice sites, the global BEC phase is unimportant. 
The advantage of using the coherent state description is 
that the Wigner function is positive. 

Due to the symmetric ordering of the expectation val- 
ues obtained from the Wigner distribution, it is difficult, 
or even impossible, to extract several correlation func- 
tions for the full multi-mode field operator, since the 
Wigner field is symmetrically ordered with respect to ev- 
ery mode. In |l3| the phase diffusion of a BEC was there- 
fore calculated by defining a 'condensate mode' operator 
associated with the projection of the stochastic field onto 
the ground state solution. Since we study the splitting of 
a BEC by a periodic optical lattice potential, it is useful 
to define analogously the ground state operators cij for 
each individual lattice site j: 

a,j(t)= dx%(x,t)ip w (x,t) , (5) 

J j th wcll 

where tpw(x,t) is the stochastic field, determined by 
Eq. 0, and i/;o(x,t) is the ground state wave function 
at time t, obtained by integrating the CPE in imaginary 
time in the potential V(x, t). The integration is over one 
lattice site. For each lattice site ground state mode cij, 
the normally ordered expectation values can be easily 
obtained (a*aj)w = + Si.j/2, etc. 

The BEC is initially assumed to be in a harmonic 
trap and we continuously increase the strength of the 
optical lattice potential until some final value, after 
which the potential is kept constant, V(x, t) = Vh(x) + 



s(t)E r sin 2 (irx/d), with s(t) — exp(nt) — 1 for t < r 
and E r = H 2 tt 2 /(2md 2 ), where d = A/2sin(0/2) is the 
lattice period, obtained by two laser beams intersecting 
at an angle 0. For very large s and close to the ground 
state only one mode per lattice site is important and the 
system can be approximated by the BHH: 

H = E ~ J (^+i + H.c.) + ^{btfb 2 } , (6) 

i 

where the summation is over the lattice sites, J ~ 
— / dxr]*(x)Crii + i(x) is the hopping amplitude between 
the nearest-neighbor sites, U ~ giD J dx\r)j(x)\ , and 
Vj = j 2 d 2 muj 2 /2, with j = site at the trap center. 
We may approximate the Wannier functions rji by the 
ground state harmonic oscillator wave function with the 
frequency lu s — 2s 1 ^ 2 E r /h at the lattice site minimum 
[2l|. When we compare the TWA results to the BHH, 
we frequently extract the expectation values involving b 
using Eq. (JSJ with b ~ a. We tested that using differ- 
ent projections does not affect the results. For n-iJ > U, 
with rii = (SjSj), the system is in the superfiuid regime 
with the long-range phase coherence and is expected to 
undergo the MI transition at m J ~ U [2^, resulting in 
a highly number squeezed ground state. 

In the numerical studies of loading the BEC into an 
optical lattice, we first solve the BEC ground state ipo by 
evolving the GPE in imaginary time in the harmonic trap 
and then diagonalize Eq. @ to obtain the quasiparticle 
mode functions Uj , Vj and energies tj . The time evolu- 
tion of the ensemble of Wigner distributed wavefunctions 
[Eq. |JIJ] is unraveled into stochastic trajectories, where 
the initial state of each realization for i/jyy is generated 
according to Eq. (J2J with the operators replaced by the 
Gaussian-distributed random variables (ay ,ot*). We in- 
tegrate Eq. using the split-step method and in several 
cases the sufficient convergence is obtained after 600 re- 
alizations. The convergence is generallyslower at higher 
temperatures. Unlike the 3D TWA 15], the ID simu- 
lations do not similarly depend on the total number of 
quasiparticle modes and we found the calculated results 
to be unchanged when we increased the number of modes. 

For the typical nonlinearity N gi D = lOOHutl, with 
I = (h/muj) 1 / 2 , the initial harmonically trapped BEC is 
well described by the GPE with the Poisson density fluc- 
tuations and the ratio between the interaction and the 
kinetic energy 7 = mgir>/(h 2 niD) < 10 -3 |2^], where 
n\D is the ID atom density. The corresponding initial 
Thomas-Fermi radius R/l = (3N g 1D /2hiul) 1 / 3 ~ 5.3. 
We take d = nl/8, resulting in E r = 32huj. Within 2R, 
we then have 30-35 lattice sites. A similar number of sites 
has also been realized in recent experiments in a cigar- 
shaped trap with d ~ 2.7^tm U In order to characterize 
the phase coherence along the lattice, we introduce the 
normalized first-order correlation function between the 
central well and its ith neighbor as Ci = |(agaj)|/ \fn^rii. 
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FIG. 1: The phase coherence between the central well and its 
nearest neighbor Ci as a function of time (left) at T — for 
different final heights of the optical lattice (curves from top 
represent s = 5, 8, 10, 15, 20, 30, 40 with the atom number in 
the central well no ~90-100). The ramping time lot — 10, ex- 
cept for s — 30, 40, lot = 6. The number fluctuations Ano for 
the central well (right) for the same runs. Here g\o = 0.05huul 
and iVo = 2000. The number squeezing can be accurately fit- 
ted according to (An ) 2 /"o ~ 0.03+0. 5e" s/8 ~ 0.03+0.3 J 4 . 



In Fig. n we show C\ and the number fluctuations 
Arii — [((a|aj) 2 ) _ (o-i^i) 2 } 1 ^ 2 m ^e central well for dif- 
ferent final heights of the periodic potential at T = 0. 
For shallow lattices the phase coherence remains high 
and steady, but for larger s it is reduced and becomes 
strongly oscillatory. Due to the large occupation num- 
bers, Ano are strongly sub-Poissonian, approaching the 
asymptotic value (Ano) 2 / no ~ 0.03 <C 1 for large s. Here 
the MI transition for the ground state is expected to oc- 
cur at s ~ 30. However, we find Ano ^ 1 for ah s, which 
can be understood by the nonadiabatic loading. 

For an adiabatic turning up of the lattice and for the 
system to remain in its ground state, we require that the 
rate of change in the tunneling amplitude to be slower 
than any characteristic time scale of the system. At low 
lattice heights it is more difficult to avoid exciting higher 
vibrational levels within one potential well, resulting in 
excitations in the higher energy bands. Moreover, the 
phonon mode energies io n in the lowest energy band de- 
crease with increasing lattice strength [3 UK and for 
high lattices it is more difficult to maintain the adiabatic- 
ity with respect to these excitations. In Fig. we find the 
number squeezing to saturate around s=20-30, indicating 
the point when an increasing number of phonon modes is 
excited and the loading becomes strongly nonadiabatic. 
Consequently, the s > 15 cases exhibit significant excess 
number fluctuations as compared to the ground state. 
After a short time period over which C± remains con- 
stant, the large An^ evolve into large phase fluctuations 
and C\ becomes oscillatory and collapses. 

The saturation of the number squeezing for strong lat- 
tices was experimentally observed in jy| for a 3D vapor in 
a ID lattice. Such a system is not tightly elongated, but 
we can still make qualitative comparisons to the experi- 
mental data. Although the saturation was assumed in 
to be an artifact of the analysis method of the interference 
measurement, we also numerically find the same satura- 
tion effect which can be explained by the nonadiabaticity 
of the loading process. If the loading is sufficiently rapid 
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FIG. 2: The coherences Ci (top left), C5 (bottom left) and 
the number fluctuations Ano (top right) for initial tempera- 
tures (curves from top) k B Ti/huj = 0, 12.5, 22.2, 33.3, 38.5 (C 5 
also with 28.5). The phase collapse time t c (bottom right) is 
evaluated at C5 = 0.5. The same system as Fig.^with s = 20. 



or the final lattice sufficiently high, so that the adiabatic- 
ity breaks down for a large number of modes, the optimal 
number squeezing is proportional to the ramping speed 
itself and the nonlinearity. Both in Fig. ^ and i n li| the 
squeezing saturates at about 15dB when riiU/J ~ 10 4 . 
The ramping time r ~ 4D00h/E r in [6j is one order of 
magnitude longer than in Fig. ^ but this is compensated 
by the weaker hopping amplitude J, so that the satura- 
tion roughly occurs at the same value of io n T. 

Although the BHH J§J is only valid for weakly excited 
high lattices, it is interesting to compare the TWA results 
to the Bogoliubov approximation to the BHH. These 
were calculated in the homogeneous lattice (yi — 0) in 
[24l I25L l2t| . Similarly, we may diagonalize the linearized 
fluctuations in Eq. JIjJ) around the ground state atom 
density with the fluctuation part 5b j — ~^2 n [f n {jd)xn — 
h^(jd)Xn], resulting in the number fluctuations in each 
site, (An^) 2 = rii Y2j l l0 j| 2 (2nj- + l) ) and the phase fluctu- 
ations between the k and I sites, (Atpki) 2 = ({<pk—<pi) 2 ) = 
V 4 Ej \rj{kd) / ^/n^ - r 3 (ld) I ^fni\ 2 (2fij + 1), where hi = 

V™iT,j( w 3Xj + w jXj) and <p t = -i/(2^/h~) Ej( r jXj ~ 
r*Xj) are the number and phase operators, with Wj = 
fj — hj, Tj = fj + hj, and fij — (x]xj)- I n the homo- 
geneous lattice with n atoms per site we have (hjj q ) 2 — 
4Jsin 2 (q(i/2)[4Jsin 2 (a d/2) + 2nU] , where q is the quasi- 
particle momentum |2J,|25j. Moreover, for nil S> J and 
N p lattice sites, (An,) 2 ~ J2 q fuj q /(2UN p )(2n q + 1) and 
(A^ fc , fc+1 ) 2 ~ J2 q fuj q /(4nJN p ){2n q + l), which at T = 
approximately yield (8n J/U) 1/2 /n and (2U/n J) 1/2 /tt, 
respectively. Numerically, we find the Bogoliubov results 
in the harmonic trap for Ano to be slightly larger and for 
At^oi smaller than the homogeneous result. The TWA 
results for An^ in Fig. are clearly larger than the ideal 
Bogoliubov limit, however, (Ano) 2 /no cx J 0A still qual- 
itatively similar to the Bogoliubov result (noli depends 
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only weakly on s). As argued in 24], if the adiabaticity 
of a phonon mode breaks down, the number fluctuations 
of the mode freeze to the value that prevails at the time 
this occurs, i.e., when Wj ~ ((t) = \d t J(t)/J(t)\. Us- 
ing the homogeneous lattice result at T — we obtain 
(An,) 2 * T, j Kj(t j )/(2UN p ). Since for all j, Q(t 3 ) is 
here roughly of the order of w, we have the asymptotic 
value for s — ► oo, (Arti) 2 ~ hiojU , qualitatively similar to 
Fig. n In order to study the effect of the nonlinearity we 
also varied in the simulations N^gm /fuol from 100 to 400 
for s = 20 and found (An ) 2 /y/nS <x U c , with c ~ —0.26, 
as compared to the Bogoliubov result c= —1/2. 

In Fig. [3 we show Ano and the coherence C\ for dif- 
ferent initial temperatures Ti for s(r) = 20. Here (Ano) 2 
increases exponentially as a function of Ti. The phase co- 
herence C5 between the central well and its 5th neighbor 
decays significantly faster than G\. The dependence of 
the phase collapse time t c on Tj is approximately linear. 
At s = 20 the effects of the harmonic trap are already sig- 
nificant, since the variation of the trapping potential over 
five sites exceeds the tunneling energy ^5 ~ 2Huj > n J. 

If the lattice potential is turned up adiabatically, the 
population of each mode remains constant and temper- 
ature T can change dramatically, as the contribution of 
each mode to T changes by the ratio of the final and ini- 
tial mode energies oj^ /lo^\ An adiabatic increase in the 
lattice strength may both increase or lower T, depend- 
ing on whether the excited band is occupied |27| and in 
the experiments the condensation temperature has been 
found to be sensitive to the lattice height pif. In Fig. [31 
we estimated the population and the 'temperature' of the 
lowest phonon modes in the TWA simulations by evalu- 
ating the projection of ipw to the Bogoliubov modes of 
the BHH © . The averages are taken over a time period 
before any significant rethermalization occurs after the 
ramping. The modes 2 and 4 are highly excited for the 
case of short t, due to the nonadiabatic loading. The ex- 
citations are damped out at higher Ti and for lot = 30, 
corresponding to W2,4T 3> 1. It is interesting to note 
that the excitations of the forth mode are only damped 
out when the rate of change in the tunneling amplitude 
£ is much smaller than the corresponding mode energy, 
or when L04 ~ 26£(t). This is more restrictive condition 
than the one found in |24j. For lot < 3, the variation 
of Ti is already completely dominated by the excitations 
due to the rapid turning up of the lattice. 

The advantage of ID lattices is that the lattice spacing 
can be easily modified by using non-parallel lasers. A 
large spacing could even allow the scattering of light, or 
the Bragg spectroscopy, from individual lattice sites and 
the separate optical detection of number fluctuations in 
each site, using a similar analysis to t 2Qj|. Moreover, an 
interference measurement on the expanding atoms can 
provide detailed information about the coherence 0. 

We studied the loading of a harmonically trapped BEC 
into an optical lattice. In a good agreement with exper- 
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FIG. 3: The contribution of the lowest modes to temperature 
(left) for ksTi/hio = 12.5 (dashed line) and 37 (solid line). 
The curves from top lot = 10, 20, 30 for both cases. The 
average temperature of the first five modes after the ramping 
lot = 30 (right). Here g 1D = 0.015hiol, N = 2000, and s = 5. 



iments |6j, we found the number squeezing to saturate 
for high lattices, which can be explained by the finite 
ramping time of the lattice potential. It is numerically 
more demanding to study a truly adiabatic loading for 
strong lattices. However, it would be particularly inter- 
esting to examine the validity of the TWA close to the 
MI ground state. Our analysis seems to indicate that, in 
the case of lattices with large filling factors, the ramping 
time required to reach the MI state may be very long and 
can be demanding in actual experiments. In the lattice 
experiments the atoms are also coupled to environment, 
resulting in dissipation with the system relaxing towards 
its ground state. We could improve our model, e.g., by 
incorporating the spontaneous emission due to the lattice 
lasers. This would introduce also a dynamical noise term 
in Eq. However, experimentally spontaneous emis- 
sion can also be avoided since, e.g., with intense CO2 
lasers the spontaneous emission rate is very low pjol ] . Fi- 
nally, our TWA studies could also be extended to the fi- 
nite temperature damping of nonequilibrium oscillations 
in a multi-well BEC what has previously been studied in 
double- well BECs H3- 
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